function [w_st,cy1_st,cy2_st]=Maxstat(gamma,betta,pi1,pi2,ey1,ey2,eo1,eo2,v1aut,v2aut)

st = @(z)[utility(ey1-z(1),gamma)+betta*(pi1*utility(eo1+z(1),gamma)+(1-pi1)*utility(eo2+z(2),gamma))-v1aut;      % IC state 1
        utility(ey2-z(2),gamma)+betta*(pi1*utility(eo1+z(1),gamma)+(1-pi1)*utility(eo2+z(2),gamma))-v2aut];       % IC state 2 

Zst = fsolve(st,[0.2 0.2]);

tau1_max = Zst(1);      % consumption of old state 1
tau2_max = Zst(2);      % consumption of old state 2
w_st = pi1*utility(eo1+tau1_max,gamma)+pi2*utility(eo2+tau2_max,gamma); 
cy1_st=ey1-tau1_max;
cy2_st=ey2-tau2_max;

end
